/*
!---------------------------------------------------------------------
! Single-center overlap functions used in Kruununhaka
!---------------------------------------------------------------------
! Copyright (C) 2010 Jussi Lehtola, University of Helsinki
!
! This program is distributed under the terms of the GNU General
! Public License as published by the Free Software Foundation;
! either version 2 of the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. In addition, the 
! author must be properly cited whenever the program or some parts 
! originated on it are used. 
!---------------------------------------------------------------------
*/

#include <math.h>
/* Hypergeometric functions */
#include <gsl/gsl_sf_hyperg.h>
/* Gamma functions */
#include <gsl/gsl_sf_gamma.h>

#include <stdio.h>

/* Maximum value for angular momentum */
#define MAXL 7
#define TABLEMAX ((MAXL)+1)
/* Table for storing normalization factors for STO-GTO overlap */
static double stogtofac[TABLEMAX];

/* Integer pow function */
double ipow(double base, int pow) {
  double result = 1.0;
  if(pow<0) {
    base=1.0/base;
    pow=-pow;
  }

  while(pow) {
    result*=base;
    pow--;
  }

  return result;
}

double powhalf(double x, int n) {
  /* Calculate x^(n/2) */
  double root;

  if(n<0) {
    x=1.0/x;
    n=-n;
  }

  if(n%2==0)
    root=1.0;
  else
    root=sqrt(x);

  return root*ipow(x,n/2);
}

/* Fill tables */
void initoverlap_() {
  int l;
  for(l=0;l<TABLEMAX;l++)
    stogtofac[l]=pow(2.0,-l/2.0-1.25)*sqrt(gsl_sf_gamma(2*l+3)/gsl_sf_gamma(l+1.5));
}

double sto_gto_overlap(double zs, double zg, int l) {
  /* Calculate overlap of normalized STO and normalized GTO */

  return stogtofac[l]*pow(zs,l+2.5)/pow(zg,l/2.0+1.25)*gsl_sf_hyperg_U(2.0+l,1.5,zs*zs/(4*zg));
}

double sto_gto_overlap_(double *zs, double *zg, int *l) {
  return sto_gto_overlap(*zs,*zg,*l);
}
